1.1 用TTree 结构存储事件数据¶
目的：学习如何将事件以TTree的格式存储到ROOT文件¶
内容：¶
- 束流入射到靶上产生核反应，靶上生成的中子和$\gamma$入射到长条型塑料闪烁体探测器的不同位置，并在不同深度与探测介质产生相互作用。
- 入射粒子形成的光信号向探测器两端传输，形成时间信号$t_u$和$t_d$, 能量信号 $Q_u, Q_d$。
- 将入射粒子的原始信息，以及探测器信息逐事件存入root文件的TTree结构中.
- 入射粒子原始信息：粒子种类、能量，粒子在探测器上的入射位置、入射深度等
- 探测器信息: 两端时间信息$t_u$、$t_d$，能量沉积信息$q_u$,$q_d$
- 打开root文件，进行基本的分析。
中子飞行时间谱测量原理¶
探测器参数¶
探测器的几何设置见下图：
- 探测器长度：2$L$, 忽略探测器宽度。
- 探测器厚度：$\Delta D$
- 光在闪烁体的传播速度: $v_{sc}$
- 光在闪烁体的衰减常数：$\lambda$
- 粒子从起始位置到探测器中心的距离：$D$
- 粒子入射位置：$x$
- 粒子相互作用的位置(深度)：$(-\frac{\Delta D}{2} , +\frac{\Delta D}{2})$
公式¶
飞行距离：
$$ d=\sqrt{D^2+x^2} $$飞行时间, 以靶为起点:
$$ \begin{equation} 中子：TOF_n(ns)=\frac{72.29824 \cdot d(m)}{\sqrt{E_n(MeV)}}\\ \gamma：TOF_{\gamma}(ns)=3.333\cdot d(m) \end{equation} $$探测器两端时间信号(不考虑由电缆长度或PMT度越时间引起的offset)：
$$ \begin{equation} t_u=TOF+(L-x)/v_{sc}\\ t_d=TOF+(L+x)/v_{sc}\\ x=(t_d-t_u)*v_{sc}/2\\ TOF=(t_u+t_d)/2 -L/v_{sc}\\ \end{equation} $$探测器两端能量，设x 处能量沉积 $2Q_0$:
$$ Q_u=Q_0 exp[-(L-x)/\lambda]\\ Q_d=Q_0 exp[-(L+x)/\lambda]\\ x=\frac{\lambda}{2}ln\frac{q_u}{q_d}\\ Q_0=e^{2L/\lambda}*Q_u*Q_d\\ $$探测器分辨率:
- 时间分辨率(FWHM)：$T_{res.}$
- 相对能量分辨率$\Delta E/E$: $E_{res.}$
随机数：
- TRandom3 class
- Uniform(Double_t x1, Double_t x2);
- Gaus(Double_t mean, Double_t sigma)
代码¶
// 运行方法：
// 将下列代码保存到tree.cc
// 在ROOT环境中运行：
// root -l //进入ROOT环境
// .x tree.cc //运行脚本
// .q //退出root环境
// root -l tree.root //打开ROOT文件
void tree(){
//常量声明
const Double_t D=500.;//cm, distance between target and the scin.(Center)
const Double_t L=100.;//cm, half length of the scin.
const Double_t dD=5.;//cm, thickness of the scin.
const Double_t TRes=1.;//ns, time resolution(FWHM) of the scintillator.
const Double_t Lambda=380.;//cm, attenuation lenght of the scin.
const Double_t QRes=0.1;//relative energy resolution(FWHM) of the scin.
const Double_t Vsc=7.5;//ns/cm, speed of light in the scin.
const Double_t En0=100;//MeV, average neutron energy
const Double_t EnRes=50.;//MeV, energy spread of neutron(FWHM)
const Double_t Eg0=1;//MeV, gamma energy
const Double_t Rg=0.3;//ratio of gamma,ratio of neutron 1-Rg
//1. 声明在tree结构中定义需要的变量分支
Double_t x;//入射位置
Double_t e;//能量
int pid; //粒子种类，n:pid=1,g:pid=0
Double_t tof, ctof;//TOF:粒子实际飞行时间，cTOF：计算得到的TOF
Double_t tu, td;
Double_t qu, qd;
Double_t tu_off=5.5;//time offset，//PMT的度越时间+电缆上的传输时间
Double_t td_off=20.4;//time offset
//2. 定义新ROOT文件，声明新的Tree
TFile *opf=new TFile("tree.root","recreate");//新文件tree.root，指针 *opf
TTree *opt=new TTree("tree","tree structure");//新tree，指针 *opt
//3. 将变量地址添加到tree结构中
//第一个参数为变量名称，第二个为上面定义的变量地址，第三个为变量的类型说明，D表示Double_t，I表示 Int_t。
opt->Branch("x", &x, "x/D");
opt->Branch("e", &e, "e/D");
opt->Branch("tof", &tof, "tof/D");
opt->Branch("ctof",&ctof,"ctof/D");
opt->Branch("pid", &pid, "pid/I");
opt->Branch("tu", &tu, "tu/D");
opt->Branch("td", &td, "td/D");
opt->Branch("qu", &qu, "qu/D");
opt->Branch("qd", &qd, "qd/D");
// histogram，ROOT文件中除了TTree结构外，还可存储histogram，graph等
TH1D *hctof=new TH1D("hctof","neutron time of flight",1000,0,100);
TRandom3 *gr=new TRandom3(0);//声明随机数
//4. 循环，计算变量的值，逐事件往tree结构添加变量值。
for(int i=0;i<100000;i++){
x=gr->Uniform(-L, L);//粒子入射位置，在-L,L范围内均匀抽样.
Double_t Dr=D+gr->Uniform(-0.5,0.5)*dD;//粒子在探测器厚度范围内均匀产生光信号
Double_t d=TMath::Sqrt(Dr*Dr+x*x);//粒子实际飞行距离，flight path
if(gr->Uniform() < Rg) { //判断为gamma入射
pid=0;
e=Eg0;
tof=3.333*(d*0.01);
}
else { //neutron
pid=1;
e=gr->Gaus(En0, EnRes/2.35); // energy of neutrons
tof=72.29824/TMath::Sqrt(e)*(d*0.01);//ns
}
tu=tof+(L-x)/Vsc+gr->Gaus(0,TRes/2.35)+tu_off;
td=tof+(L+x)/Vsc+gr->Gaus(0,TRes/2.35)+td_off;
ctof=(tu+td)/2.;//simplified calculation.
hctof->Fill(ctof);
//neutron：energy of recoil proton in plas. q0=0-En； gamma：q0=0-Egamma，compton plateau
Double_t q0=e*gr->Uniform();
qu=q0*TMath::Exp(-(L-x)/Lambda);
qu=gr->Gaus(qu,qu*QRes/2.35);//QRes, relative energy resolution
qd=q0*TMath::Exp(-(L+x)/Lambda);
qd=gr->Gaus(qd,qd*QRes/2.35);
opt->Fill();//5.将计算好的变量值填到Tree中
if(i%1000==0) cout<<".";
}
cout<<endl;
// 6.将数据写入root文件中
hctof->Write();//写入预定义的histogram到文件
opt->Write();//写入TTree到文件
opf->Close();//关闭文件
}
练习：¶
- 理解上述1.-6.步骤的逻辑关系。自行练习如何将数据已TTree格式存到ROOT文件中。
- 理解如何在模拟中加入探测器分辨率信息。
查看ROOT文件¶
- juypter运行tree->Draw()时，需要有网络链接才能正常显示。
In [1]:
//%jsroot on
In [2]:
TFile *ipf=new TFile("tree.root");//root -l tree.root
ipf->ls()//ROOT 环境下 > .ls
In [3]:
//观察tree的结构
tree->Print()
观察指定事件数据¶
In [4]:
tree->Show(0);//显示指定的事件数据 event=1
In [5]:
tree->Show(100);
In [6]:
TCanvas *c1=new TCanvas();//* 在ROOT环境下可省略
c1->Clear();//* 在ROOT环境下可省略
tree->Draw("td-tu>>htx(500,-20,50)");//位置一维分布
c1->Draw();//* 在ROOT环境下可省略
In [7]:
tree->Draw("tof>>hh(1000,0,100)");//实际飞行时间
c1->SetLogy();
c1->Draw();
In [8]:
TH1D *hh=(TH1D*) ipf->Get("hctof"); //在代码中得到文件内数据指针的方法，在ROOT环境下课直接用 hctof->Draw()
hh->Draw();//显示文件内histogram
c1->Draw();
思考题1：观察上两个图有什么不同，分析其原因¶
In [9]:
c1->SetLogy(0);
gStyle->SetPalette(1);
tree->Draw("ctof:td-tu>>(2000,-20,50,1000,0,100)","","colz");//二维图显示
c1->Draw();
检查参数之间的关系是否与物理条件一致¶
- gamma的到达探测器的时间(飞行时间)与探测器的x位置无关？
In [10]:
tree->Draw("ctof:td-tu>>(2000,-20,50,50,40,45)","pid==0","colz");//加入条件,gamma
c1->Draw();//为什么在两侧上弯？
In [11]:
tree->Draw("td-tu:x","","colz");//观察两个参量的关联
c1->Draw();
In [12]:
tree->Draw("log(qu/qd):x","","colz");
c1->Draw();
In [13]:
!jupyter nbconvert 1.1_create_tree.ipynb --to html
In [ ]: